Reformulation of the LDA+t/ method for a local orbital basis 
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We present a new approach to the evaluation of the on-site repulsion energy U for use in the 
LDA+U method of Anisimov and collaborators. Our objectives are to make the method more 
firmly based, to concentrate primarily on ground state properties rather than spectra, and to test 
the method in cases where only modest changes in orbital occupations are expected, as well as for 
highly correlated materials. Because of these objectives, we employ a differential definition of U. 
We also define a matrix U, which we find is very dependent on the environment of the atom in 
question. The formulation is applied to evaluate U for transition metal monoxides from VO to NiO 
using a local orbital basis set. The resulting values of U are typically only 40-65% as large as values 
currently in use. We evaluate the U matrix for the e g and ti g subshells in paramagnetic FeO, and 
illustrate the very different charge response of the e g and tig states. The sensitivity of the method 
to the choice of the d orbitals, and to the basis set in general, is discussed. 
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I. INTRODUCTION 

The understanding and evaluation of the electronic 
structure of strongly correlated materials is a long- 
standing problem. For weakly correlated materials such 
as nearly-free-electron-like metals, covalent semiconduc- 
tors, ionic solids, and even rather complex intermetallic 
transition metal compounds, the local density approxi- 
mation (LDA, which we understand to include the spin 
degree of freedom as well) to the exchange-correlation 
functional that occurs in density functional theory gives 
very reasonable ground state properties and even band 
structures (which are excited state features). For cor- 
related materials, however, LDA can be completely 
wrong: the now-classic example is the canonical cuprate 
La2CuQ4, which LDA predicts to be a non-magnetic 
metaoH whereas it is actually an antiferromagnetic insu- 
lator. Model many-body Hamiltonian treatments, such 
as the Hubbard model,0 can readily explain the observed 
type of ground state but do so in terms of adjustable 
parameters and the neglect of many aspects of the crys- 
tal that may influence most of its properties. Evaluation 
of the dynamic self-energy, which gives the description 
of excitations, is appropriate for comparing with many 
experiments, but even low-order approximations can be 
very tedious to evaluated 

Within the past few years Anisimov and collaborators 
have proposed an extension of the LDA approach (now 
called LDA+ti) based on lessons learned from Hubbard 
model studiesta that single out a particular local orbital 
and the associated on-site repulsive interaction U as the 
fundamental characteristic to be addressedErQ. They pro- 
pose that the LDA treats the effects of U reasonably well 
in some average sense, even in highly correlated systems, 
but that one must allow a deviation from this average 
behavior by including a correction to the total energy 
including a term like 
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where J is the exchange constant and n ms is the new 
charge that includes a local charge redistribution (rel- 
ative to the LDA value n rns ) obtained by solving the 
LDA+t/ equations self-consistently. The local orbital 
and spin indices are m and s, respectively. It is assumed 
in the method that one can identify the orbitals to be 
treated (d orbitals of Cu in La2Cu04 for the example 
mentioned above). 

The LDA+C7 method achieves some spectacular suc- 
cesses, such as leading to an antiferromagnetic insulating 
state of La2CuC>4 with band gap and atomic moment 
in reasonable correspondence with observed valuesu, and 
leading to similarly impressive descriptions of the transi- 
tion metal monoxides. There remain questions, however, 
such as the proper way to specify the orbitals, the cor- 
rect way to obtain the interaction constants (U and J), 
and how, if possible, to extend the method to give an 
improved treatment of the metallic phase when the in- 
sulator is heavily doped. In this paper we address these 
questions. A primary feature is that, since the method 
is perforce focussed on an atomic orbital, it is natural to 
use a local orbital basis set. We will refer to the local 
orbital of interest as the "d orbital(s)" although in some 
applications they may be / or, rarely, s or p orbitals. 



II. DESCRIPTION OF LDA+U AS CURRENTLY 
PRACTICED 

In extending the LDA method to account for corre- 
lations resulting from strong on-site interactions, there 
are several criteria that one might hope to satisfy, such 
as (1) it should reduce to LDA when LDA is known to 
be good; (2) the energy is given by a functional of the 
density; (3) the method specifies how to obtain the local 
orbital in question (perhaps through a self-consistency 
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procedure); (4) the definition of U and J are provided 
unambiguously; (5) the method predicts antiferromag- 
netic insulators when appropriate; (6) the description of 
highly correlated metals is improved over the LDA de- 
scription. This list, although perhaps still incomplete, is 
already very ambitious, and only certain of these desires 
have been addressed seriously. _ 

Anisimov, Zaanen, and Andersen (AZA)Q chose to re- 
fine the LDA by including an orbital-dependent one- 
electron potential to explicitly account for the impor- 
tant Coulomb repulsions not treated fully in LDA. This 
was accomplished in analogy with Hartree-Fock theory 
by correcting the mean field contribution of the d — d 
on-site interaction with an intra-atomic correction. This 
correction has been applied in slightly varying forms, but 
a representative example of the functional to be solved is 



Elda+u = E LDA [n] - ^2Ni[n ma ](Ni[n ma ] - 1) 

i 
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Here i denotes the lattice site, and terms involving J 
have been neglected because we do not need to specify the 
complete form of the functional for this paper. iVj is the 
site sum of the d charges, evaluated for the self-consistent 
LDA+C/ densities. The second term is presumed to be a 
reasonable description of the direct Coulomb interaction 
energy included in the LDA expression. 

Equation (2) reveals that, in the LDA+J7 approach, 
one singles out beforehand the atomic orbitals to be 
treated, and decides how to specify them. Implementa- 
tions to date use orbitals arising in the linearized muffin- 
tin orbital (LMTO) method. The d orbitals to which 
the U correction is applied are numerical solutions to a 
Schrodinger equation inside an atomic sphere, and are 
zero outside this sphere. In addition, the LDA+J7 is 
clearly no longer a straightforward density functional be- 
cause it depends on parameters U and J that depend on 
the LDA density rather than the LDA+C7 density. 

The one-electron potential is the conventional LDA 
form of potential, plus an orbital-dependent shift of en- 
ergy given by 



Av ms = ua 







(3) 



if Umm' — * U is orbital independent. The changes in the 
electronic structure are proportional to U, and the defi- 
nition and calculation of U is the next topic to address. 

To obtain U and J, AZA performed LMTO calcula- 
tions for a supcrcell in which the d charge on-one atom is 
constrained and the eigenvalue is obtained. E3 The d or- 
bitals on all atoms in the supercell are decoupled entirely 
from the remaining part of the basis set. This makes the 
treatment of the local orbitals an "atomic-like" problem, 
which greatly reduces the difficulty associated with con- 
straining the occupation numbers. It also has the effect 



of leaving a rather artificial system to perform the screen- 
ing. For example, in NiO the screening system consists of 
oxygen p orbitals that cannot hybridize with the Ni d or- 
bitals, plus whatever other virtual orbitals are included 
in the basis set. 

The discreteness of the d eigenvalues makes it simple 
to specify the charge in the spin-orbitals in the supercell, 
and U and J are determined from the relations 



U - £ 3 <iT (f + 2J f) ~ £ 3dT (I + 2' § _ !) 



(4) 

in which the d occupation differs by unity around a mean 
polarization of unity, and 

J = £3dT (f + 2"' § - 3) - £ 3di (f + §, f - 3) , (5) 

which is a straightforward difference between up and 
down eigenvalues for unit spin polarization. Here 
£3^1- (n-t, fit ) ( £3dj,( n T' n l)) ^ s t ne spin-up (spin-down) 3d 
eigenvalue for occupancies nj and n^. 

While it is widely recognized that the on-site repul- 
sion U is a screened quantity, the manner in which the 
screening should be done is not precisely -specified. An 
early study by Cox, Coulthard, and LloycO for 3d met- 
als used a renormalized neutral atom approach, although 
it was recognized that screening processes might extend 
over a somewhat larger region. Anisimov and collabo- 
rators have chosen the method presented in this section, 
but in this paper we pursue a different approach, with a 
somewhat different objective. 



III. REFORMULATION OF LDA+U FOR A 
LOCAL ORBITAL BASIS 

We specify in following subsections the various ways in 
which our approach differs from that in current use. 



A. LCAO Basis Set 

We begin with a basis set of local orbitals {(/)}, whose 
lattice sums lead to the standard linear combination of 
atomic orbitals (LCAO) Bloch basis functions for the 
one-electron Hamiltonian.E3 To represent an occupied 
atomic orbital (including core states) we use a contracted 
set of Gaussian functions, multiplied by appropriate an- 
gular functions for s, p, or d behavior. In particular, we 
choose at the beginning (from a neutral atom or an ion) 
the d orbitals of central interest. Although we have no 
indication of any better choice than the d orbital of the 
corresponding atom (e.g. neutral Cu in La2Cu04), our 
method allows the ability to check how sensitive the re- 
sults are to the form chosen for the orbital. In addition to 
basis functions describing filled atomic (or ionic) orbitals, 
we add other Gaussian functions to the basis to provide 
a more nearly complete basis for the valence and conduc- 
tion states than a minimal basis set would provide. This 
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feature is an advantage of our local orbital representa- 
tion, as the ability to include self-consistent screening by 
a crystalline density of general form in the calculations 
is important. 

This LCAO basis set brings up an important feature. 
As a sum of squares of wave functions, the charge den- 
sity contains two types of terms. One consists of atom- 
centered contributions containing the coordinate depen- 
dence <t>t m (r — R) and is clearly identified as a contribu- 
tion to the charge from angular momentum I of the atom 
located at R. The other contribution has a coordinate 
dependence given by 4>i m (r - R)(j>i> m >{r - R'), R ± R', 
which at a particular point may be positive or negative 
and cannot be assigned uniquely to any atom. The Mul- 
liken decomposition^ which assigns half of each such 
term to charge component I, m on the atom at R and 
the other half to charge I'm' of the atom at R', is widely 
used when atomic decomposition of the charge is desired. 
Mullikcn population is well understood to be not only ar- 
bitrary, but also dependent on the flexibility of the basis 
set, and therefore should not be endowed with any im- 
portant physical meaning. 

A central fact that must be addressed is that the total 
charge density cannot be decomposed, precisely or mean- 
ingfully, into simple atomic contributions alone. This fact 
means that the orbital occupations Tii^rn 

that are the cen- 
terpiece of the LDA+i7 approach unfortunately are not 
particularly well defined. For our LCAO basis set we will 
use for n ms charge contributions solely of the first type, 
which will be called on-site charges to distinguish them 
from Mulliken charges. These on-site quantities also can- 
not properly be called occupation numbers since there is 
no sum rule for their total, and it is not impossible that 
for a given orbital the value can exceed unity. 

B. Specification of the functional 

Although we do not carry out LDA+U calculations 
in this paper, we are thinking in terms of a generalized 
LDA+U functional that is consistent with our philoso- 
phy behind the correction. Without more formal justi- 
fication than is normally done in the LDA+J7 approach 
(and which we do not address seriously here) , any change 
must simply be tested to see if it produces better results. 
The form that we envision has affected our study of how 
to define and to evaluate the interaction constants that 
arise in the method. We suppose that the correction is 
to provide adjustment to full potential LDA results, and 
therefore includes both a suborbital index and a spin in- 
dex on the reference charges n — > n ms . These numbers 
will differ, sometimes greatly, for different irreducible rep- 
resentations of the point group of the atom. The correc- 
tion then might be written suggestively as 

EhDA+U = -ElDa[«] + ^ (Um,m' ~ S s , s 'J m ,m') 



This change may affect the types of orbitally-ordered so- 
lutions that will be obtained. This form ensures that 
the LDA solution is an exact stationary solution of the 
LDA+U functional (for which the correction vanishes 
identically), which is not the case for Eq. (2), i.e. if 
shell-averaged values of n are used. Aside from strongly 
correlated solids, another interest of ours is to ascertain 
whether LDA+U can provide a useful improvement of 
the description of 'simpler' systems such as the transi- 
tion metals Fe and V, where anisotropy (relative amounts 
of tig and e g character) is not reproduced accurately in 
LDA, or in correlated metals where no bandgap occurs 
but charge rearrangement might be appropriate. 

C. Procedure for determining U and J 

We take as our ansatz that the constant U (resp. J) 
occurring in the LDA+U functional should describe the 
cost in potential energy of charge (resp. spin) fluctua- 
tions in the actual crystal, i.e. with all normal interac- 
tions and degrees of freedom available to the electrons. 
Thus we do not decouple the d states from surrounding 
states. The U and J terms are going to be applied in pre- 
cisely the same system from which they are determined. 
We comment below on the question of dealing with the 
associated cost in kinetic energy due to charge fluctua- 
tions. For the remainder of this paper we concentrate 
solely on U, postponing a related treatment of J such as 
suggested by Solovyev et alt3 for the future. 

We also take the point of view that the main purpose 
of LDA+C/ theory, as in density functional theory, is to 
obtain ground state properties, rather than to approx- 
imate excitations with the eigenvalues. Describing the 
ground state may require small rearrangements of occu- 
pation numbers away from their LDA values, and usually 
less than one-half, so we employ a differential definition 
of U (also used by Solovyev, Dederichs, and Anisimovtj) 
rather than one employing occupation numbers differing 
by unity. We will see that this introduces extra richness 
into the charge rearrangements described by the LDA+Z7 
method, because a small change in (say) ti g population 
can be strongly compensated by a change in e g popula- 
tion. 

We employ then a generalized constrained depsity 
functional approach as proposed by Dederichs et. aJl3 to 
calculate the change in energy due to constraints on lo- 
cal orbital densities. We minimize the local density func- 
tional subject to the constraint that on-site local orbital 
charges n Q S be equal to designated values Q a<s , where 
a labels an irreducible representation of the point group 
of the atom in question (e.g. t2 g or e g ), and that total 
charge N be conserved: 

£{Q)= min {£; LDA [n T , ni] 
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fi{ / n(r)d 3 r- N)}. 



(7) 



in energy due to constraining a set of orbital densities in 
the manner of Eq. (7). Since there is no change in energy 
if the charges are "constrained" to be their LDA values 
Q = Q LDA (so q = 0), the energy change is given by 



The Lagrange multipliers are the usual chemical poten- 
tial /i and the potential shifts w a , s necessary to satisfy 
the constraint n a , s — Q a . s . Dependence on the total 
number TV of electrons (always conserved) will not be 
displayed explicitly. Variation with respect to the spin- 
orbitals leads to a one-electron Schrodinger equation in 
which the potential is the LDA potential, supplemented 
by local orbital shifts w ajS on the orbitals in the irre- 
ducible representation a having spin s. These additional 
shifts of potential can be represented as a non-local po- 
tential 
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(8) 



a.s m£a 



where {<fi} are normalized atomic orbitals. 

Evaluation of the constrained energy in Eq. (7) de- 
serves comment. Solution for the constrained energy in- 
volves generating the Kohn-Sham equations, which have 
an additional potential of the form of Eq. (8) that effec- 
tively constrains the density as desired. The conventional 
method of evaluating the energy is to sum the resulting 
eigenvalues and correct for double counting of the Hartree 
energy and the miscounting of exchange-correlation en- 
ergy. That cannot be done directly, because the Kohn- 
Sham eigenvalues contain the effects of the additional 
potential of Eq. (8) and one does not obtain Elda- The 
additional term that has been included by summing the 
LDA+[/ eigenvalues however contains only the additional 
one-body term ^ Q w aiS n atS , and this term can be sub- 
tracted to obtain -Elda evaluated for the constrained den- 
sity. 



D. The Constrained Energy 



dq 



(10) 



subject only to the condition that £^ is analytic (as we 
assume) . 

The general behavior of the constrained energy can 
be seen by noting that w is linear for small changes in 
occupation, i.e. linear in q. Since at the minimum of Eq. 
(7) n = Q, we may use these quantities interchangeably 
to write 



w 



where Sf 



-XJq + O(q) 2 = -VSn + 0{Sn) 



(11) 



?LDA 



and U is the constant (matrix) of 
proportionality. For the remainder of this section we con- 
cern ourselves with the linear "response" that is implicit 
in the LDA+/7 method, although we demonstrate in the 
numerical results of Sec. V.A where non-linear correc- 
tions begin to arise. Then the energy shift is given by 



where 



u 



l 



-q- U • q, 



-dw/dq. 



(12) 



(13) 



The constrained energy £{q) can be decomposed into 
the kinetic energy term, the interaction with the external 
potential, and the remainder, the potential energy: 



% = 



£q,K + £q,cxt + £q,P- 



(14) 



£<3,cxt is linear in q and gives no contribution to U, but 
the quadratic term involving U contains both a kinetic 
energy contribution Uk and a potential energy contribu- 
tion Up, 



It is convenient to introduce a vector notation for the 
local occupations, the constraining values, and the (La- 



grange parameter) potential shifts: 



and sim- 



ilarly for Q and for w. Since we will be dealing with 
quantities relative to their LDA values, we also use the 
notational conveniences 



q=Q-Q 



LDA 



£ q - = £(Q)-S(Q LDA ). 



From the Hellman-Feynman-likc relation 



dq 



Vq£q 



0) 



we can generalize the constrained-density functional the- 
ory viewpoint of Dederichs et al£B to obtain the change 



U 



u. 



u 



K(P) — ^7q~Vq£q,K(P)- 



(15) 



In a self-consistent calculation, any change in local or- 
bital charge results in an accompanying change in ki- 
netic energy as well as a potential energy change. Thus 
one might argue that it should be Up that goes into 
the LDA+U calculation, and the kinetic energy change 
in the constrained LDA calculation should be removed: 
U — > Up = U - U K is the appropriate "t/" in LBA+U. 
Apparently this is the correction (in our language) that 
Anisimov and GunnarssonO expected to account for by 
disconnecting their local orbitals from all other basis 
functions. 

Using the Hcllman-Feynman relation Eq.(9) to obtain 
w(q) it is straightforward to obtain U of Eq.(15). The 
kinetic energy contribution Uk can be evaluated directly, 
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which is a very delicate task numerically, or alternatively 
from the relation £k = — [d£ / d\ogm] m=mo which makes 
use of the variational nature of the total energy £ (m D 
is the electron mass). Results for Uk will be presented 
elsewhere. For now, the value of U that we evaluate and 
report below is the total value U=Uk+Up. 



E. Change in Independent Variable 

It will be instructive to consider the potential shifts w 
to be the independent variables in an associated energy 
functional leading to q(w) rather than w(q). This is also 
in keeping with the practice in the constrained density 
approach of choosing the shifts w and then calculating 
the charge response q. This change of variable is done by 
a Legendre transformation 



u-0 



£w = £q + q- w, 

which from the differential forms 

8Eq = — w ■ Sq 
=> 5£yj = q ■ Sw, 

leads to the energy shift 

£$ — / q(w) ■ dw 
Jo 
1 



(16) 



(17) 



-w ■ U 1 • w. 



(18) 



This formalism brings in the matrix U _1 implicit in 
Eq. (11) relating the charge shifts in various suborbitals 
to potential shifts applied to other suborbitals, e.g. a de- 
composition of the Hubbard U for d orbitals into e g and 
tig contributions for cubic site symmetry. This result is 
reminescent of the extension of of the definitions of U and. 
J [Eqs. (4)-(5)] by Solovyev, Hamada, and Terakurata 
to give different values U e and Ut 2g , but their procedure 
did not provide off-diagonal terms. The effects of differ- 
ing charge response in the e g and ti g channels will be 
quantified in Sec. V. The concept can be extended to 
non-site-diagonal interactions, viz. d orbitals interacting 
with neighboring oxygen p orbitals. 

We now establish a sum rule relating the matrix ele- 
ments of U to the conventional scalar U, which for clar- 
ity we denote Udd = dwd/dQd, where Qd is the total d 
charge and Wd is a shift in potential applied to all d or- 
bitals. Since a change in potential Wt 2 acting on the t^ g 
orbitals followed by a change in potential w e acting on 
only the e g orbitals is equivalent to a potential Wd of the 
same magnitude acting on all d orbitals, we have, in the 
linear regime 







d 







dw t2g dw 
By definition rid = nt 2g 



dw d 

n e , so from the definition 



(19) 



dn a 



(20) 



we have a sum rule relating the matrix elements to the 
conventional Coulomb repulsion constant 



a,/3=t 2g ,e g 



(21) 



We provide below a numerical test of this sum rule for 
NiO. 



IV. METHOD OF CALCULATION 

For the metallic constituents of the compounds we con- 
sidered, a basis set representing six s-, four p-, and three 
rf-type functions is expanded on a set of sixteen Gaussian 
functions. The O basis set is expanded on a set of twelve 
Gaussian exponents contracted into four s- and three p- 
type functions. The Coulomb and exchange-correlation 
potentials comprise the effective potential, V e fr, which 
is also described by a superposition of atom-centered 
Gaussian-type functions. By choosing this expansion, the 
matrix elements of the Hamiltonian are analytic. Details 
of the method, and comparison to results of the full po- 
tential linearized augmen±fid_.plane wave method, have 



tential linearized augmented-.]: 
been published elsewhere £3'Ej 



For this work it is important to obtain sufficiently well 
converged values of orbital densities. Tests using special 
point meshes in the irreducible 1 /48 of the simple cubic 
Brillouin zone (IBZ) for eight atom cells up to 56 k points 
indicated that ten or twenty /c-points in the IBZ gave the 
necessary accuracy. A temperature broadening of 0.07 
eV was used to facilitate convergence to self-consistency, 
and it was verified that this size of broadening did not 
change the results. 



V. EVALUATION FOR TRANSITION METAL 
MONOXIDES 

We have applied this approach to evaluate U for the 
transitions metal monoxides MO, M=V, Mn, Fe, Co, 
and Ni, in the paramagnetic state and for the cubic rock- 
salt structure. The (experimental) lattice constants used 
were: VO, 4.093 A; MnO, 4.444 A; FeO, 4.332 A; CoO, 
4.260 A; NiO, 4.193 A. 



A. Sub-orbital Independent U 

First, applying a potential shift Wd equally to all d sub- 
orbitals analogously to LMTO treatments, the derived 
value of U is shown in Table I. Comparison is provided 
with values obtained by the method of AZA, and it is 
seen that the values obtained are 40-65% of the values 
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TABLE I. Calculated values of U for transition metal ox- 
ides, compared to values of Anisimov, Zaanen, and Andersen 
(AZA) .13 Empirical values include representative values from 
the literature. 



Ref. VO MnO 


FeO 


CoO 


NiO 


This work 2.7 3.6 


4.6 


5.0 


5.1 


AZA 6.7 6.9 


6.8 


7.8 


8.0 


Empirical 4.0-4.8° 7.8-8.8° 


3.5-5.1° 


4.9-5.3° 


6.1-6.7° 


7.0 d 


3.9 c ,7.0 d 


4.9 C 


7.9°, 6.1 c , 7.5 d 



a V. I. Anisimov et at, Ref. 6. 

b M. R. Norman and A. J. Freeman, Phys. Rev. B33, 
8896 (1986). 

C J. Zaanen and G. A. Sawatzky, J. Solid State Chem. 88, 
8 (1990). 

d A. E. Bocquet et al, Phys. Rev. B46, 3771 (1992). 

TABLE II. d shell charges, according to various definitions, 
for transition metal monoxides from VO to NiO. 



Type 


VO 


MnO 


FeO 


CoO 


NiO 


Formal 


3 


5 


6 


7 


8 


On-site 


3.67 


5.45 


6.22 


7.41 


8.22 


Mulliken 


3.09 


5.48 


6.44 


7.20 


8.40 



obtained by AZA. That our values are smaller is no sur- 
prise, since our approach (of not disconnecting d orbitals 
from other orbitals) naturally allows additional screening 
to occur, by including hybridization between d orbitals 
and neighboring oxygen p orbitals. Moreover, charge re- 
arrangement between e g and ti g subshells reveals that 
there is some intra-d-shell screening in the current ap- 
proach. In addition, our definition of the d orbital is not 
identical with that of AZA. 

Figs. (l)-(3) for MnO, FeO, and NiO illustrate the 
change in subshell charge with Wd as well as the total 
change, which is what determines U. Taking MnO, for 
example, it is seen in Fig. 1 that the effect of positive 
Wd is to decrease n t2g as expected, but n Sg instead in- 
creases. Clearly charge rearrangement within the d shell 
is leading to a reduction in U (i.e. additional screen- 
ing). Similar behavior occurs for FeO (Fig. (2)) while 
for NiO the e g charge remains almost unchanged as Wd 
is varied. Note that, besides the differences in approach 
and in basis sets, the values obtained in the AZA ap- 
proach are evaluated for differences in d charge of unity. 
Empirically determined values (obtained by comparing 
to excited state data) lie in nearly all cases between our 
values and those of AZA. 

The values of U in Table I are obtained as the first 
derivative of polynomial fits to the Wd vs. Qd curve 
[Eq. (pi])]. Figs. (l)-(3) indicate the AQ^ vs. Wd curve 
for shifts w d up to ± 0.8 eV for MnO, FeO, and NiO. 
The change in total d charge is linear up to this size of 
shift («20-30% of U). Even for this size shift, however, 
the individual ti g and e g contributions are beginning to 
become nonlinear, as seen most clearly for MnO in Fig. 1. 

On-site charges and Mulliken charges within LDA, for 



0.5 



O 

u 
o 

<1 



0.3 



0.1- 



-0.1 



-0.3 



-0.5 



X 


MnO 

„+ 


A"' 


Sum 








\^ 












X 



1.0 -0.5 



0.0 0.5 

w d (eV) 



1.0 



FIG. 1. Change in on-site d charge (solid line) in MnO 
resulting from a potential shift applied to all d states on 
a single Mn atom in a four molecule supercell. The total 
charge is decomposed into its e g (short dashed line) and t2 g 
(long dashed line) components. 
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FIG. 2. Change in d charge in FeO, plotted as in Fig. 1. 
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our basis set, are compared in Table II. The charges are 
less ionic than their formal (dipositivc) charge, as expe- 
rience would suggest. (Although atomic charge within a 
crystal cannot be defined uniquely, it is widely accepted 
that 'effective' ionic charges are nearly always reduced 
by hybridization from their formal, full ionic values.) Al- 
though VO is somewhat of an exception, the Mulliken 
charge does not differ more than 4% from the on-site 
charge for these examples. The response of Mulliken and 
on-site charges are very different, however, with Mulliken 
charges varying more slowly. If one uses Mulliken charges 
rather than on-site charges to obtain U, the resulting val- 
ues are much larger: 3.8 eV for VO, 6.2 eV for MnO, and 
11.1 eV for NiO. 



B. Sub-orbital Dependent U 



a 
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We have studied FeO charge redistribution when e 
and t2 g subshclls are treated separately. In Fig 
present the change in subshell charge when a shift in 
potential is applied individually to the subshclls. In both 
cases, charge forced out of one subshell by an upward 
shift in potential goes primarily into the other subshell, 
amounting to very strong intra-d-shell screening in these 
cases. Using Eq. (20), we obtain, in eV -1 , 



r-l 

t 2 g J 2 g 
-1 



U e g ]t 



1.18 
-1.39 



-1.00 
1.41 



(22) 



which satisfies the sum rule of Eq. (21). The inverse is, 
in eV, 
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FIG. 4. Change in subshell (e g and tig) charge in FeO, 
resulting from a potential shift of only one of the subshells. 
Top panel: change in tig charge; bottom panel: change in e g 
charge. The label indicates the type of applied potential shift 
w: d indicates shifts of all d states; tig (e g ) indicates shift of 
only tig (e g ) states. Only positive energy shifts are shown. 
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Recall that C/ rf(i =4.6 eV, so in the usual orbital- 
independent treatment the corresponding matrix would 
be 

(U t2g , t2g U t2g , e \ _ ( 4.6 4.6 \ 

Thus the behavior that looks rather peculiar in Fig. 4, 
and the negative off-diagonal elements in Eq.(22), do not 
lead to pathological behavior in the direct matrix U. 

C. Dependence of Local Orbital Shape 

The LDA+C/ procedure is built around some choice of 
local orbital. In an LCAO basis, this orbital is specified 
at the beginning, and we have used neutral atom d or- 
bitals from atomic LDA calculations. Another possible 
choice might be, say, the d orbital from a positive ion. For 
FeO we have checked the effect of using the Fe 2+ d orbital 
obtained from an atomic calculation on an isolated ion. 
The difference in radial density is shown in Fig. 5. For 
the FeO (paramagnetic) solid, the on-site charge of 6.22 
electrons (Table II) changes to the rather peculiar value 
of 4.85 electrons, and the calculated value of U increases 
from 4.6 eV to 7.8 eV. The total energy, however, changes 
only by +0.12 eV/FeO, which is a very modest change 
(adding / functions in an LCAO or LMTO calculation 
can result even larger changes, which are unimportant for 
most purposes) . It is clear that the choice of d orbital can 
affect the calculated value of U, certainly in the LCAO 
method but most likely in any calculational approach. 



VI. DISCUSSION 

The results presented in the previous section reflect a 
strong difference in response of the e g and ti g electrons, 
at least to potentials of moderate strength. Such dif- 
ferences have been noted several times in the literature. 
In the context of the LDA+Z7 method, Solovyev et al.i3 
have advocated using using separate values of U for the 
two subshells in perovskite structure transition metal ox- 
ides. Their method of obtaining U a was a generalization 
of the standard method described in Sec. II. We, on the 
other hand, have adopted the differential definition of U 
that leads to a matrix U a p. 

To begin to understand the response of the separate 
subshells, we show in Fig. 6 the on-site e g and ti g den- 
sities of states (DOS) on a Ni atom in an eight atom 
supercell of NiO, both before and after a downward shift 
of all d states by -0.54 eV. The t2 g states in the rock- 
salt structure are weakly dpn bonding and form a narrow 
band, whereas the e g states form dpa bonds that produce 
wider e g bands. From Fig. 3 it is seen that such shifts 
produce negligible change in the on-site e g charge, with 
all the difference coming entirely from the on-site ti g sub- 
shell. This result is counterintuitive, since the ti g DOS 
is full, and pulling it down seemingly cannot increase its 
occupation. The e g DOS is open shell and could accept 
charge, but does not do so. 

The resolution of this paradox lies in the change in 
the representation of charge of the system by the LCAO 
basis functions as a shift in potential is applied. By look- 
ing at other local DOS, for both on-site and Mulliken 
charge decompositions, we have found that a downward 
shift of d states, which changes the degree and character 
of hybridization as well as the probability of occupation, 
results in a more active participation of the virtual or- 
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bitals in representing the charge density. To some extent 
this is a real effect: d — p hybridization increases as the d 
states are pulled down nearer the p bands, and what one 
calls the d function, or the d charge, becomes less well 
defined. (The definition becomes clear for well separated 
atoms, and perhaps for states well separated in energy 
from any other states.) The effect is present with other 
basis sets as well, but is more difficult to identify. 

It is important to understand clearly the source of this 
paradox. It arises because the "d charge density" is not 
a precisely and objectively defined quantity. Although 
we are accustomed to thinking in terms of five d bands 
in a transition metal oxide that can be identified and 
whose DOS integrates to five electrons per spin, this is a 
fiction that becomes apparent as soon as orbital overlap 
becomes appreciable. This difficulty has the same origin 
as the difficulty in defining the "d orbital" to be used 
in the LDA+t/ method. These ambiguities are problems 
that must be lived with until a better prescription can 
be formulated. 



VII. SUMMARY 

We have presented a reformulation of the method of 
obtaining U for an LDA+U calculation. The approach 
is based on a local orbital expansion, which is a natural 
one considering that the d orbital is to be singled out and 
specified anyway. We aim specifically to improve ground 
state properties rather than to account for spectroscopic 
data. 

Values of U using this approach are found to be only 
40-65% of the values of Anisimov et al. Most of this 
difference is understood in terms of the definitions and 
procedures that are used in each case. A novel feature 
here is the identification of an interaction matrix that 
describes interactions that are non-diagonal in the sub- 
orbital index, e.g. the change in energy of ti g states due 
to a change in e g charge. The off-diagonal parts of this 
interaction are expected to be strongly dependent on the 
environment of the ion, and this expectation is borne out 
in our study of FeO. 

There are important aspects of our approach that re- 
quire further work. The contribution to U from the ki- 
netic energy, and how it should be dealt with, is one loose 
end. The most appropriate choice of d orbital is another 
question that may require some experience to answer. 
Carrying out LDA+[7 studies to compare with results 
using the previous LDA+f/ method, and ascertaining the 
effect of off-diagonal interactions, are however the main 
priority, and this work is in progress. 
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